;************************************************
; NCL Graphics: 
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"

;************************************************
begin


  load "load_setup.ncl"

  vname  = "ABSORB"
  vunits = "" 
  vscale = 1. ;;e-3

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  fa = ca+ma
  fb = cb+mb
  casa  = (/"A: "+na, "B: "+nb/)
  caseM  = (/pa+fa,pb+fb/)
  case  = (/caseM(0),caseM(1)/)
  Ncase = dimsizes(case)
  plotname = "figure_"+vname+"_zonal_"+shortname

  print(" ")
  print(" ")
  print("..................................................................")
  print(" ")
  print(" ")
  print(" varname   : "+plotname)
  print(" ctrl      : "+case(0))
  print(" test      : "+case(1))
  print(" plot name : "+plotname)
  print(" ")
  print(" ")

;************************************************
; Open netCDF file and read in dimensions
;************************************************

   filename=case(0)+".nc"
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)
   ;print("nlev "+sprinti("%7.2i",nlev))
   ;print("nlat "+sprinti("%7.2i",nlat))
   ;print("nlon "+sprinti("%7.2i",nlon))
   
;************************************************
; set default res
; res1 for gsn_csm_contour
; resp for gsn_panel
;************************************************
  wks   = gsn_open_wks (plotfmt,plotname)             ; open workstation
; gsn_define_colormap(wks,"gui_default")          ; choose colormap
  gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")          ; choose colormap
  res1                      = True                ; plot mods desired
  res1@gsnFrame             = False               ; advance page
  res1@gsnDraw              = False                
  res1@gsnSpreadColors      = True                 ; use full colormap  
  res1@gsnSpreadColorStart  = 30                    ; start at color 10
  res1@gsnSpreadColorEnd    = 235                    ; end at color 96
  res1@cnFillOn             = True                 ; turn on color
  res1@cnLinesOn            = False                ; turn on/off contour lines
  ;res1@cnFillMode           = "AreaFill"           ; default "AreaFill" "RasterFill"
  ;res1@cnFillMode           = "RasterFill"           ; default "AreaFill" "RasterFill"
  res1@lbLabelBarOn          = True               ; turn on/off the labelbar
  res1@lbOrientation         = "Vertical"        ; labelbar Orientation "Vertical"
  res1@lbLabelsOn            = True               ; turn on/off the label of labelbar
  res1@lbBoxMinorExtentF     = 0.40                ; area allotted to each box of the LabelBar    
  res1@lbLabelFontHeightF    = 0.02                                               
  res1@lbLabelAutoStride     = False                ; optimal labels         
  res1@cnLineLabelsOn        = False                ; turn on line labels
  res1@cnLevels             = (/1.e-8,2.e-8,5.e-8,1.e-7,2.e-7,5.e-7,1.e-6,2.e-6,5.e-6,1.e-5/) 
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  res1@cnLevelFlags         = False
  res1@tiXAxisFontHeightF    = 0.020               ; XAxis Font
  res1@tiYAxisFontHeightF    = 0.018               ; YAxis Font
  res1@tiYAxisString         = "Approximate Pressure(hPa)"
  res1@tiYAxisString         = "Pressure(hPa)"
  ;res1@tiXAxisString         = "Latitude"
  res1@vpXF                  = 0.12     ; x location left
  res1@vpYF                  = 0.88     ; y location  up
  res1@vpWidthF              = 0.65     ; width
  res1@vpHeightF             = 0.35     ; height
  res1@gsnStringFontHeightF  = 0.02
  res1@tiMainString          = " "  
  res1@tiMainFontHeightF     = 0.03
  res1@tiMainFontThicknessF  = 0.35
  res1@trYReverse            = True
  res1@gsnYAxisIrregular2Linear = True

  resP                 = True          
  resP@gsnPanelYWhiteSpacePercent = 1  ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 5  ; draw panel with white space added
  resP@txString            = "" 
  resP@gsnPanelLabelBar    = False                ; add common colorbar
  resP@lbLabelFontHeightF  = 0.012               ; labels size
  resP@lbBoxMinorExtentF   = 0.30                ; area allotted to each box of the LabelBar 


  res1@tmYRMode             = "Automatic"          ; turn off special labels on right axis

;************************************************
; read in variables & create plot
;************************************************

  plot = new(4,graphic)
  
  do i=0,Ncase-1
     filename=case(i)+".nc"
     fi=addfile(filename,"r")
    ;print(" read "+case(i))

     freqi = fi->FREQI 
     var   = fi->$vname$ 
     PS    = fi->PS

     freqi = where(freqi.le.0.,-999,freqi) 

     fvar  = var  
 
     P0mb =1000.
     hyam = fi->hyam
     hybm = fi->hybm

     ; type of interpolation: 1 = linear, 2 = log, 3 = loglog
     interp = 2
     ; is extrapolation desired if data is outside the range of PS
     extrap = False

     ; create an array of desired pressure levels:
    ;pnew = (/ 1000.,925.,850.,800.,750.,700.,600.,550.,500.,450.,400.,350.,300.,250.,200.,150.,100./)

     ;; use cam5 l30 levels, find and use the standard levels 

     pnew = (/ 100., 121., 143.,  \ 
               168., 200., 233., 274., \ 
               300., 379., 446., 500., \ 
               610., 700., 763., 821., \ 
               850., 887., 913., 936., \ 
               957., 976., 993., 1000. /) 

     pnew = pnew(::-1) 

     ; note, the 7th argument is not used, and so is set to 1.
     pvar = vinth2p(fvar,hyam,hybm,pnew,PS,interp,P0mb,1,extrap)

    ;printVarSummary(pvar) 

     var3d = dim_avg_Wrap(pvar(time|:,lev_p|:, lat|:, lon|:)) ;;dim_avg_Wrap(fi->$vname$(time|:,lev|:, lat|:, lon|:))

     var3d!0 = "time"
     var3d!1 = "lev_p"
     var3d&lev_p = pnew
     var3d&lev_p@units = "hPa"
     var3d!2 = "latitude"
     
     var2d = var3d(0,:,:) 
     var2d@_FillValue = -999

    ;printVarSummary(var2d) 
    ;print(var2d&lev_p@units) 

     var2d=var2d*vscale
     res1@gsnLeftString = casa(i)
     res1@gsnRightString= vunits                


    ;print(var2d) 

     var2d = where(var2d.le.0.,-999,var2d) 


    ;plot(i) = gsn_csm_contour(wks,var2d(:,:), res1) 
     plot(i) = gsn_csm_pres_hgt(wks,var2d(:,:), res1) 

     if(i.eq.0) then
        vara = var2d 
     end if
     if(i.eq.1) then
        varb =  var2d
     end if

     delete(var2d)
     delete(var3d) 
     delete(hyam) 
     delete(hybm) 
     delete(freqi) 
     delete(var) 

  end do 



  P0mb =1000.
  hyam = fi->hyam
  hybm = fi->hybm

  vara = where(vara.le.0.,-999,vara) 

  varc = vara 
  vard = vara 
  varc = 100.*(varb - vara)/vara
  vard = varb - vara


  res1@gsnSpreadColors      = True                 ; use full colormap 
  res1@gsnSpreadColorStart  = 20                    ; start at color 10
  res1@gsnSpreadColorEnd    = 195                    ; end at color 96


  res1@gsnLeftString = "(B-A)/A"
  res1@gsnRightString= "%"
  delete(res1@cnLevels) 
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  res1@cnLevels             = (/-50,-30,-20,-10,-5,0,5,10,20,30,50/) 
 ;plot(2) = gsn_csm_contour(wks,varc(:,:), res1) 
  plot(2) = gsn_csm_pres_hgt(wks,varc(:,:), res1) 


  res1@gsnLeftString = "B-A"
  res1@gsnRightString= vunits 
  delete(res1@cnLevels)  
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  res1@cnLevels             = (/-5.e-7,-2.e-7,-1.e-7,-5.e-8,-2.e-8,-1.e-8,0.,1.e-8,2.e-8,5.e-8,1.e-7,2.e-7,5.e-7/) 
  plot(3) = gsn_csm_pres_hgt(wks,vard(:,:), res1) 





  resP                 = True
  resP@gsnPanelYWhiteSpacePercent = 1  ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 1  ; draw panel with white space added
  resP@txString            = ""
  resP@gsnPanelLabelBar    = False                ; add common colorbar
  resP@lbLabelFontHeightF  = 0.012               ; labels size
  resP@lbBoxMinorExtentF   = 0.30                ; area allotted to each box of the LabelBar 


  
  gsn_panel(wks,plot,(/2,2/),resP)
  
end

